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Abstract 

We investigate the effect of the superlattice structure on the single particle 
transport along the c-axis of high temperature superconductors. In partic- 
ular, superlattice systems that consists of metals/insulators and d-wave su- 
perconductors (NS/IS superlattice) are considered. We find that for the NS 
superlattice in the large mass anisotropy limit of the metal, the density of 
state in the low energy section is bulk d-wave like except that the position of 
the quasi-particle peak can be considerably smaller than the gap value, while 
for the IS superlattice, the quasi-particle peak remains at the gap value. We 
also calculate the nonlinear differential conductance in the planar junction 
measurement. It is found that the width of the Andreev peak at zero-bias 
may be affected strongly by the superlattice structure, specifically, it can be 
considerably reduced due to the destructive interference of the Andreev re- 
flections from all the superconductors. Such a reduction in the width makes 

the Andreev peak resonant-like and has been observed in a recent experiment. 
PACS numbers: 74.20.-z, 74.50.+r, 74.80.-g, 74.25.Fy 
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I. INTRODUCTION 



The artificial prepared superlattices have been an important system for investigating the 
relevancy of periodic structure to the physical properties of the bulk system. The idea is that 
if the dephasing length of the electron extends over many periods, the superlattice structure 
becomes relevant. Therefore, by engineering the periodic structure appropriately, one is able 
to tune the effective parameters that electrons experience, thus it may enable us to access 
more regions in the parameter space. In particular, one may be able to engineer the density of 
state so that the transport properties can be changed. Such systems have motivated a lot of 
studies both experimentally and theoretically 0. There are many outstanding issues arising 
from the investigation of superlattice systems. For instance, the giant magnetoresistance 
(MR) observed in magnetic superlattices, such as Fe/Cr, has triggered serious investigation 
of the transport through heterostructure of metal and ferromagnetic materials . 

The superlattice that has superconductors involved is another important system. It was 
shown in an early study that certain superlattice structure can close the energy gap of s- 
wave superconductors |J, which is interesting by itself and may have real applications in 
the future. The interest in supperlatices with superconductors has been further boosted 
up since the discovery of the high temperature superconductors (HTSC) [|]. The layered 
structure along the c-axis and particularly the existence of closely packed CuO^ planes makes 
HTSC a natural-occurring superconducting superlattice system. A great body of work has 
thus been devoted to investigate the c-axis dynamic properties during the past 0. After 
the intense work for so many years, it now becomes clear that in the overdoped region, 
the normal state of high Tc superconductors is metallic like for both c-axis and in-plane 
directions ||. In particular, the splitting due to hopping between adjacent Cu0 2 planes 
in a unit cell are seen directly in recent photoemission experiments ]7|. While these works 
imply that single particle hopping definitely occurs along the c direction, it is not clear how 
exactly one should model the transport properties of the high Tc cuprates along the c-axis 
in the superconducting phase and how the superlattice structure affects the single particle 
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transportation along the c direction. 

The splitting observed in the photoemission experiments implies that neighboring Cu0 2 
planes in a unit cell are strongly coupled by a large hopping integral. Therefore, one can 
model it by a single superconducting layer, which is coupled to the superconducting layers in 
adjacent unit cells via an effective hopping integral t e ff. This is also consistent with the fact 
that the coherence length along the c direction is about the width of each superconducting 
layer. The problem is then how one should model layers between adjacent superconducting 
layers. It was pointed in Ref. || that there exist two distinct limits, depending on the ratio 
t e ff to T c . In the limit when t e ff is much smaller than T c , the system behaves more like 
a NSNSN ...NSN (or a ISISISISI) superlattice, i.e., a normal metal-superconductor (NS) or 
insulator-superconductor (IS) superlattice. On the other hand, if t e ff is close to T c in mag- 
nitude, different unit cells are coupled strongly and the system behaves as a usual anisotropic 
superconductor. In this work, we shall be mainly investigating the nonlinear transport in the 
limit when t e ff is much smaller than T c . In order to model the c-axis transport in HTSC, 
we shall consider several concrete models, for instance, the NS and IS superlattice models 
as both of them can give rise metallic behaviors in c and ab directions in the normal state. 
Another possible modification is to consider a NS superlattice with a large mass anisotropy 
(effective mass in the c-axis << effective mass in ab directions) in the metal so that elec- 
trons essentially hop along the c-axis in the NS superlattice model. Two quantities will be 
calculated: (1) the density of state for single particles, which will be useful for the transport 
measurement with high resistance contact. (2) the differential conductance (dl/dV) for 
measurements with low resistance contact , such as measurement made by planar junctions. 
We shall include the d-wave nature and consider general Fermi surface topology. Our results 
will be useful for artificial superlattices such as the YBa 2 Cu 3 7 ^s/ Pr Ba 2 Cu 3 7 ^s super- 
lattice system Q]. At phenomenological level, they also provide insights to the issue that to 
what extent the c-axis of HTCS can be thought as a NS superlattice or a IS superlattice, in 
particular, how the superlattice structure along the c-axis affects the tunneling spectrum. 
Our results indicate that in the case for a NS superlattice, if the mass anisotropy is small, 
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the superlattice structure will induce subgap structure in the density of state. For large mass 
anisotropy in the metal, however, the bulk d-wave-like density of state is reproduced with 
the quasi-particle peak shifted to a smaller value, in contrast to the case when we model 
the C~ clXIS clS du IS superlattice where the quasi-particle peak remains at the gap value. The 
NS superlattice model provides a possible explanation of why some measurements of the 
gap size along the c direction give smaller values @||. We also calculate the dl/dV curve 
for measurements using planar junction or point contact with large contact area along the 
c-direction. Two important features are found to be due to the superlattice structure: (1) 
The width of the Andreev peak at zero bias is considerably smaller than the gap due to 
destructive interference of Andreev reflections from all the superconductors. This gives a 
natural explanation to why in some old measurement, the Andreev peak along the c-axis 
appeared to be so sharp and was thought as zero-bias resonant conductance peak |1(J. (2) 
There exist distinct oscillations in the region when V is larger than the gap size. This is 
also a result of interference from all in the interfaces. In particular, the IS superlattice will 
provide a rising background so that these oscillations gradually rises as V increases. Both 
features are observed in a recent measurement of Au/Bi2212 junctions near T c \TL . 



This paper is organized as follows: In Sec. II, the density of state is calculated in the 
framework of Bogolubov de-Genne equation. We show that in the large mass anisotropy 
limit of the NS superlattice, the usual d-wave like behavior is reproduced. We also calculate 
the case when it is a SI superlattice and briefly compare these results. In Sec. Ill, the 
non-linear dl/dV curve for planar junctions is calculated and its relation to experiments is 
discussed. Finally, we give a concluding remark in Sec. IV. 

II. THEORETICAL FORMULATION AND DENSITY OF STATE 

We shall start by investigating the density of state for a NS superlattice. A NS superlat- 
tice with a spherical Fermi surface and s-wave for superconductors was previously considered 
by Hahn M to model transport along the c-axis for YBCO. To account for the c-axis trans- 
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port of high Tc cuprates, we shall extend it to include d-wave nature and use a more general 
Fermi surface. In particular, during our calculation, we shall also briefly address the effect 
due to a particular form for the c-axis hopping t c = —t±cos 2 (2<f)). This particular form of 
hopping is suggested from the measurement of angle-resolved photoemission (ARPES) and 
band theoretical calculation [0,|T^]. Here is the angle tan~ 1 (q y /q x ) (see Fig.l). 

We shall first calculate the density of state. In the most general case, the density of state 
is given by 



n{E)= j J jdj0dn 5{E _ E{qx ^ K)l 



(1) 



when k is the Block wavenumber along the c-axis. For given E, q x , and q y , if K(E,q x ,q y ) 
can be found, n(E) can be rewritten as 



n(E) = JJ 



dq x dq y , „> 



9re 



(2) 

. In this case, 



where we have defined the one dimensional density of state by n\{E) 
one needs to find K(E,q x ,q y ). For a NS superlattice, a suitable framework for calculating 
k(E, q x , q y ) is the Bogolubov de-Genne equation which can be written as 



(3) 



where h = £ q — fip, A = Ao cos(20) in the superconducting part while in the metal part, we 
simply set Aq = 0. We shall assume a general form for £ 9 



h 


A 




u(r) 


= E 


u(r) 


A* 


-h 




v(r) 




v(r) 



h 2 h 2 

= (q 2 + q 2 ) H q 2 in the superconductor, 

2m' y 2ms 



2^(?x + Qy) + 7^-Qz in the metal. 



2m N 



Since the system is translationally invariant in x and y directions, q x and q y are con- 



served so that we can write u(r) 



mz 



i(q x x+q y y) 



Eqs.(|^) then reduces to one dimensional equations. In the metal part, A = 0, q x 
and q y are on the Fermi surface so that we may parameterize q x = kp sin 9 cos (f> and 
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q y = kp sin 6 sin 0, where kp = \J 2 ™z F and 9 is the azimuthal angle along the c direc 
tion. In this parameterization, for a given energy E, the particle (+) and hole (— ) momen 
turn become 
jfe ± = 



— sin 2 9) ± \JE 2 — A 2 ] in the superconducting part, while 



2mjy 



(//p cos 2 6 ± E) in the metal. Furthermore, Eq.(f2|) becomes 



n(£) 



kp sin 6 1 cos 9d<pd9 8k 



(4) 



(2tt) 3 a^' 

It is important to note that the factor 1 — m/m' sin 2 9 has to be positive as it represents 
(p^) 2 at the Fermi surface in the absence of A. Therefore, large m/m', i.e., large mass 
anisotropy in the metal, restricts electrons to hop only along the c-axis. 

The wavefunctions u(z) and v(z) have to to continuous at boundaries. Since the system 
is periodic, it is sufficient to impose the boundary conditions within a unit cell (see Fig. 1) 

»jf(0)=iis(0),%(0)=i; s (0) 
u' N (0)=u' s (0),v' N (0)=v' s (0) 
u s (b) = Xu N (-a),v s (b) = Xv N (-a) 

u' s (b) = Xu' N (-a),v' s (b) = Xv' N (-a) (5) 

where the last two equations follow from the Bloch theorem with | A |= 1. Following Ref. 
13fl and using the above boundary conditions, one finds that A satisfies 



D + DxX + D 2 X + D 3 X + D A X 







(6) 



where D , Di, .., _D 4 are all real. Since| A |= 1, the four roots to Eq.(||) are in the form 
eX p±(«0i) an( j eX p±( i6, 2)_ This implies that Do = D^,D\ = D3. Thus the solutions to Eq.(^|) 
is completely determined by the ratio and 



cos(k d) 



Do J D 



(7) 



where we have expressed A by exp(mcf) with d = a + b. We shall see that the positive 
sign is for the particle excitation, while the negative sign is for the hole excitation. It is 
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convenient to measure the length by £jv = h 2 k F /y/mm]srA an d energy by A (e = E/A ). 
Denoting A = a/^, B = b/£ N and e F = fi F /A , after some algebra and using the Andreev 
approximation |13|], we find that for a fixed <fi and 6 



_A = -4 [cos(2e F AX) cos(2e F7j BX m )r 1 - sin(2e F AX) sin(2e F 7 J BX m )r 2 ] , (8) 
Do 

^ = 4T 2 2 + (g 2 - 1) [1 + cos(4e F AX) + cos(4e F7j BX m ) - 2r 3 ] 

+ (g 2 + 1) cos(Ae F AX) cos(4e F 'yBX m ) - 2gsm(4e F AX) sm(4e F -fBX m ). (9) 



Here X = cosO, 7 = ^Jms/rriN-, X m = \Jl — (m/m 1 ) sin 2 9 and 5 = (cos6*/7X m + 
7X m / cos#)/2. For e 2 > cos 2 20, Tj are given by 

Fl = cos (t) cos (x^) - ^ sin (t) sin ( £ iB ' (10) 

r2 = ^ cos (t) cos (xf ) - ^ sin (t) sin © ' 

T 3 = cos 2 [1 + cos(4e F AX)] + cos 2 [1 + cos(4e F T 5X m )] , (12) 

where = y/e' 1 — cos 2 20; while for e 2 < cos 2 20, one simply replaces cos(sin ) by cosh 



(sinh), and by y/ cos 2 20 — £ 2 . 

Homogenous case This is the case when m = m' and m^v = %. Therefore, g> = 1. 
We have Ti = T 2 = T(e), as a result, we find that cos(kgT) = rcos[2e F (A + B)X] ± 
Vl — T 2 sin[2e F (74 + B)X], recovering results obtained in Ref. |L3| with k+ ~t k and p+ ~t p 



being approximated by 2e F (A + B)X. It is instructive to check the case of A = 0. In 
this case, we have p = and thus T = 1. We find that k = k . Therefore, ± represents 
particle and hole channels respectively. For the homogeneous case, since T 2 < 1, we may 
write T = sin a(e). It is then easy to see that n\{e) = ^ ±^ are the same for particles and 
holes. As shown in Fig. 2(a), if we use negative e to represent the hole, the one dimensional 
density of state is an even function of e. Fig. 2(b) shows the result when ms 7^ mjv, it is 
seen that there is no particle-hole symmetry any more. Note that the values of A and B are 
chosen to be close the data for YBCO: £n ~ £,ab ~ 1.5nm, a ~ 0.85nm, and b ~ 0.38nm. 
It is also instructive to check that the NS superlattice model also includes the limit when 



t e ff is close to T c in magnitude as one can simply set A = 0, B = 1 in the homogeneous 
case. In this case, T is simply cosh ^Je 1 — cos 2 20 for e 2 > cos 2 20 and cos Vcos 2 20 — e 2 for 
£ 2 < cos 2 20. Thus rii(e) = \ , E = for e 2 > cos 2 20, which reproduces the usual d-wave 

Y e 2 — cos 2 2</i 

density of state after integration over 0, see Fig. 3 (a). 

Fig. 3(b) shows a case when A becomes nonzero. It is obvious that the introduction 
of any small N section (i.e., small A) moves the quasi-particle peak into the subgap region 
while leaving in a peak-like structure at E = Aq. Fig. 3(c) shows a more realistic case for 
A = 0.57 and B = 0.21 (homogeneous). The structure at E = A becomes too small to be 
observed. 

Large mass anisotropy limit We now discuss the effect of increasing the ratio m/m'. A 
possible candidate to describe the high Tc cuprates is the limit when m/m' becomes so large 
that electrons essentially can only hop along the c-axis between CuC>2 planes in different 
unit cells. As we pointed it out already that large m/m' restricts 9 to be small. For infinite 
m/m', we can simply set X = cos 9 = 1 and perform the following reduced integral 



n r (E) 



kpdcf) dn 



(13) 

. Consider the case when m$ = 



x=i 



(2vr)3 dE x , 

If the superconducting part is s-wave, n r (E) is p^rff 
m^, it is easy to see that Xg = g = 1. We find that Ti = T 2 = r(e) and cos(k ± 2ef(A + 
B)X)d = r(e). The point is that the one dimensional density of state dn/de starts from a 
gap whose size is determined by the root to r(e ) = 1 for Eq < 1. Near Eq, one has 



l-l = 1 1 \ dT \ « 1 J T '^ 1 (14) 
'de' dv/T^T 2 '^ d\j 2 v/F^V v ; 

which has a BCS-like square root singularity but at smaller energy (eq < 1). Therefore, in 
the lower energy section, we obtain a bulk s-wave like density of state with quasi-particle 
peak at smaller energy. For high energy sections, near the zone boundary of k, the BCS-like 
square root singularity will repeat again. For d-wave superconductors, Eq.(|T3"D is replaced 
by 



dK 1 1 8T / P(£ ) I COS 20| 

l de l dv/T^p'cfe'^V 2 ^-^08 201' 1 J 



where Sq depends on 0. Thus, the BCS-like square root singularity is preserved for a fixed 
0. In Fig. 4, we show the result of direct integration over for different ms/m^ ratios. It is 
clear that at low energy section, n r (E) behaves similar to the bulk d-wave density of state, 
except that the quasi-particle peak moves to smaller energy. If we take Ao ~ 30 — 40 meV, 
the quasi-particle peak is around 8 — 11 meV which is close to what experiments have seen 
|§]. Note that another difference between the above results and the true bulk d-wave density 
of state may lie at high energy section. In general, the superlattice structure introduces 
coupling among the wavevector of the electron and the reciprocal lattice vectors of the 
superlattice, i.e., (0, 0, 2wk/(£) (= Q). The coupling strength depends on the magnitude 
of the Fourier component A (Q). If the superlattice periodicity is good, A (Q) will not 
be small. The superlattice structure will force n r (E) to have similar d-wave like structure 
whenever the zone boundary (cos kg? = ±1) is encountered in high energy sectors. However, 
if the system is not large enough along the c direction or the dephasing length of the electron 
is short, one may expect that A (Q) is small and hence the repeated structure in high energy 
section will not appear. 

We now address the effect due to the anisotropy of the c-axis hopping. In particular, 
we consider t c = — t± cos 2 (20). As it is clear from the above calculations and also the 
following calculation for the IS superlattice, the density of state depends only on the ratio 
of ins t° m N (defined as 7 2 ). Obviously, if both the intra-superconducting cells and inter- 
superconducting cells c-axis hoppings follow the same form. There will be no effect at all. 
However, new features may arise if they do not follow the same form. In particular, we find 
that if either m s or m N does not depend on 0, a second peak could arise inside the subgap 
region. In views of past experimental findings, this seems to be unlikely. Therefore, we shall 
not consider such possibility in the following. 

IS superlattice As a comparison, we now consider the case when the materials between 
Cu02 planes in different unit cells are modeled by insulators. For this purpose, we introduce 
a large potential V to every metal cell so that the metal cells effectively become insulators. 
We shall consider the simplest case when m = m! = mjy = %. When the superconducting 

9 



gap vanishes, this model reduces to the standard Kronig-Penny, therefore, it represents 
a natural generalization of NIN superlattice. In this case, becomes purely imaginary. 
Following the same procedure, we find that Eqs. (H) and (pf) are essentially the same except 



that g is now defined by i(a + l/a)/2 with a = y /xf cos 2 9/ (V — /xp cos 2 9) and cos(2£pAX) 
and sm(2eFAX) being replaced by cosh(2epyiX/a) and sinh(2ei?y4X/a). The new Tj are 
given by (for e 2 > cos 2 20) 



r a = cosh (^) cos (^) + ig^ sinh (^) sin (^) , (16) 
T 2 = -ig cosh (^) cos (^) - L sinh (^) sin (^) , (17) 



cos 2 



1 + cosh(- 



a 



2 /eAa 



+ cosh 2 ( — ) [1 + cos(4£ F BX)) , (18) 



X 

while again for e 2 < cos 2 20, one replaces cos(sin ) by cosh (sinh), and by y/ cos 2 20 — e 2 . 
In Fig. 5, we show a typical result for the three dimensional density of state for the SI 
superlattice. Because for energy below the gap, e 2 < cos 2 20, quasi particles and quasi 
holes are both evanescent, the quasi-particle peaks now move to the gap value (except that 
there is a slight asymmetry in particles and holes). Technically, this is due to that for 
e 2 < cos 2 20, all factors in Tj are not oscillatory and thus the right hand side in Eq.(0) is 
greater than one, without any propagating solutions. This is similar to the case for a bulk 
d-wave superconductor as we analyzed for Fig.3(a), so one gets quasi-particle peaks right at 
the gap value. 

The dramatic difference in the positions of the quasi-particle peaks for the NS and IS 
superlattices offers a useful check on both models. For artificial made NS superlattices, the 
reduction in the position of the quasi-particle peak has been observed in Ref. H. However, 
for the c-axis measurements on real high Tc materials, the results are controversial and 
remains to be clarified in the future |T3|. 
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III. NONLINEAR DIFFERENTIAL CONDUCTANCE 



We now investigate the differential conductance (dl/dV) for the measurement with low 
resistance contact, such as measurements made by planar junctions or by point contact with 
large contact area. The schematic setup is shown in Fig. 6. In principle, at first and last 
interfaces that connect with the electrodes, there may exist insulating layers. This will be 
included as interfacial scattering potentials. The most important feature coming out from 
this type of measurement, in comparison to the high contact resistance measurement, is the 
possible exhibition of the excess current due to the Andreev reflection (AR). However, for a 
long time, the AR along the c direction was not observed until quite recently Andreev peaks 
with reduced widths are observed in Au/Bi2212 junctions near T c ||1 1|| . 

In this section, we shall carry out a theoretical calculation of the dl/dV curves based on 
the NS/IS superlattice models. We shall show that the width of the Andreev peak in such 
systems is often reduced due to the superlattice structure. A similar problem in which the 
superconductor in each cell is replaced by another metal (N') was analyzed in Ref. [15[ using 
a full quantum mechanic approach. Here we shall approximate the voltages in each layer 
(either in N or S) by constants such that the slope of the voltage, thus the electric field, is 
fixed. This is exact in the superconducting cell, however, in the metal cell the approximation 
is valid only if the electric field is weak (specifically, eVa/nd has to be smaller than ep). 
To obtain the differential conductance, we need to match the quasi-particle wave functions 
at boundaries: z = nd and z = nd + b. At the boundary z = nd, the quasi-particle wave 
function for the normal metal is 

u(z = nd) \ [ 1 ^ 

ip(z = nd) = = A r , 



u(z = nd) 
v(z = nd) ^ 



Ak„nD 



(A 



B r , 



( \ 



ik„nD 



e ik "" D + D, 



V 1 / 



V / 
V 1 / 



-ik„ nD 



(19) 



where we have measured energy in unit of A and lengths by hence D = d/^ and 
fc± = J cos 2 6 + (±£- V*)/e F with V* = ^-(1 - n/N ) and N being the total number of 
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layers. For the superconducting region at the right hand side of z = nd, the wavefunction is 



Here 



[z = nd) = E n+1 



+G. 



n+l 



(u \ 
yv OJ 



e ipt +1 nD + p 



n+l 



the 



f u ' 
yv OJ 
U \ 



\ U °J 



e -ipt +1 nD 



(20) 



wavevec- 



tor p n+ i is measured by kp and is given by yji(l — 72 sin 9) + / yi/ep(±e^ — V^+i)] with 



7i 



and 7 2 = V^ +1 = ^[1 — (ri-D + B)/N D] (measured by hp), and the coherent 



factors Uq = (1 + e^j£)j2 and 



i-'n 



^1 — s ( j ) /e)/2. The boundary conditions at z = nd are 



ip(nd) = (p(nd), 



mN 9 a( a\ 9 it a\ 2m NH . , 

W^ 0M ~ d- z ^ nd) = -j?~* (nd) ' 



(21) 



where we have introduced a delta potential with strength H at this interface |TJ|. This 
provides a relation between A n ,B n , C n , D n and £7 n +i, F n +i, G n +i, H n +i, which is conveniently 
expressed by introducing a transfer matrix M\ 



Mi 



( F ^ 

Fn+l 



Here the full expression of Mi is given by 



( A ^ 

-Tin 



Br, 



Cm 



y -D n y 



2 ^ fen 7l fc n 2 ^ fen 71 fc n 

mo h _ j 2Z — p l+i] e i2e F[pt+i+ k i} nD — i 2Z + p l+^] e - i2e F[pt +1 - k i] nD 

2 ^ fe? 71 2 71 &n 

2 L fc~ 7l fc n 2 fen 7l fc n 

Bl[X _ — ?l+L] e i2e F[Pn +1 +kn]nD voh— l 2 ^. _|_ p "+i 1 c ~ i2€ F fp^ +1 -fen ]w-P 



(22) 
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2 ' feJF TlfcrP 2 ' fe»F Tl^rF-' 



1 J2e F [p~ +1 +k+]nD t^M _ , P n +i l c ~i2e F \p^ +1 ~k^]nD 

2 ' fe« 71 fen J 



Bin _ ;M _ p n+i 

2 [i V 7lfcn- 

]nD «aM _|_ j2Z_ _ Pn± 

2 ' 7ifc. 



2 



Pn+l j c -»2e F fp„ + 1 +fc n 1nD 



f[l + ip + %je . 

21 fen 71 fen J 21 fen 71 fen J 

_ Pn+l 1 c i2e F [p~ +] +fc~1nJ> «oh ^2Z _|_ Pn+l l c -»2e F [p~ +1 -fc'jnD 

fen 71 fen'' 2 L fen 71 fc n / 



(23) 



with Z (= fsfp-) characterizing the interface potential. Note that in principle different 
layers may have different Z. The above formalism is for the NS superlattice. To model a IS 
superlattice, one adds a large potential V to the metal cell so that k^ become evanescent 
and are replaced by ^Jcos 2 9 + ^(±e — — ^)- 

Similar boundary conditions are also imposed at z = nd + b, and yield the relation 



M 2 



1 F ^ 



r n+l 
Gn+1 
H n+ \ 



1 A ^ 



(24) 



where M 2 can be obtained by Mi simply by the following changes: nD — > nD+B, k^ — > 
and if — > —if. Thus the effective transfer matrix that connects the nth cell and (n+l)th 
cell has the form 



M, 



( A ^ 

B n 

C n 

\ Dn J 



( , \ 

B n +i 

D, 



'n+l 



(25) 



where M n ^ n+ i = M 2 ■ Mf 1 . Note that coherence in superconducting region is assumed so 
that no phase memory is lost when quasi-particles tunnel from one superconducting cell to 
another one. The transfer matrix of the whole system for N - layers, T No , is thus obtained 
by 



iVo-l 
n=0 



(26) 
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so the coefficients in the N -th cell A^ , B^ , Cn ,Dn are connected with those in first cell 
A , B , Cq, D by the relation 



T 



N 



B 
Co 
K D 0J 



^ An ^ 
Bn 
C No 

\ D N j 



(27) 



We shall further impose the boundary conditions near electrodes. Firstly, because electron 
can only flow into the system from the high voltage electrode, we require 



4) = i,A) = o. 



Secondly, in the low voltage electrode, we impose the condition 



(28) 



B 



N 



c 



No 



0. 



(29) 



so that no particle (both hole and electron) flows into the system from the lowest voltage 
side. Eqs. ( p8|) and(^) are then sufficient to determine the amplitude B and C in terms 
components of T No , we find that 



r} / rp31 rp23 rp21 rp33 \ I ( rp21 rp33 rp32 rp23 ^ 

n — \ 1 N 1 N ~~ 1 No 1 N ) I \ 1 N 1 Nq ~ 1 N 1 N ^ 



S < / rp3\ rp22 rp2\rp32\ I / rp23 rp32 rp22 rp33 N 

°0 — \ 1 N 1 N ~ 1 No 1 No) I \ No No — 1 N 1 No^ 



(30) 



To obtain the differential conductance in terms of the above coefficients, we apply the same 



argument used in Ref. flnj in which one treats the superlattice between electrodes as a 
scattering center. We obtain the same expression for the differential conductance 



dl/dV = 2N{E F )ev F A{\- \ B (eV) | 2 + | C (eV) 



(31) 



where A is the area for the junction. 

To illustrate the behavior of the differential conductance, it is useful to consider the sim- 
plest case when Nq = 1, i.e., a NSN structure. This will be the more realistic configuration 



for measurement on what is often referred as the NS junction in the literature [16|. We shall 
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consider the case when m = m' = = ms and ep » 1 so that one can apply the so-called 
Andreev approximation in which one approximates and p^ by 1 (6 = and n = 0), i.e., 
the wavevectors along c-axis are approximated by hp. This is the approximation used to 
derive the differential conductance of the NS junction in the BTK theory |I6|| . Nevertheless, 



it turns out that for a NSN junction, such approximation leads to a result with zero Andreev 
current. A closer examination shows that in the large Sp, one has to keep the expansion of 
both and p^ to the first order in (±£ — V^)/ep or (±£^ — V^ +l )/ep when computing the 
phase terms. We find that such manipulation leads to 

B = 

— I _|_ e ^4> Bi 

C = 2~ 2 2e,Bi U V 0- ( 32 ) 

-Uq + VQe 2e t- B% 

Therefore, for e 2 < cos 2 20, is purely imaginary and hence in the limit of large B, Co 
approaches — vq/uo, recovering the standard BTK result. For e 2 > cos 2 20, is positive 
and C oscillates with periods, determined only by B. The oscillation of C was actually 
observed in Ref . . It is a result of interference of Andreev reflection from both interfaces 
(left and right interfaces) of the superconductors. In Fig.7(a), we show the dl/dV curves of 
directional tunneling (0 = and without using the Andreev approximation) for different B 
with A being fixed at one. One sees that the period of oscillation decreases for increasing 
B, in consistent with the above analysis. 

As iVo exceeds one, the lengthscale A also participates in determining the oscillations. 
In this case, the Andreev approximation always yields Bo = when Z = because Tj^ 
and T^ 2 are always zero in this approximation. The total conductance is thus determined 
by Co given by —Tj^/T^. Fig. 7(b) shows the calculated differential conductance for the 
directional tunneling (electrons incident perpendicular to the interface) of a NS superlattice 
by using the exact expressions of the transfer matrix. As a demonstration, here we consider 
the homogeneous case with the total number of layers being ten. One sees that a repeated 
main Andreev peak is already seen around e ~ 4, resembling the repeated structure when 
crossing the zone boundary in infinite systems. Note that these repeated structures may 
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not be exhibited in real experiments as the dephasing length could easily get shorter in 
higher voltage. Nevertheless the main Andreev peak at zero-bias will still be a special 
feature to superconducting superlattices. Most importantly, unlike the width of the Andreev 
peak in the NS junction, the width of the Andreev peak can be much less than Aq. This 
is due to destructive interference of the Andreev reflections from all the superconductors. 
As the normal metal diminishes, the separation between these Andreev peaks increases, 
and eventually, only the one at the zero bias survives. At the same time, the width of the 
Andreev peak extend to Ao- In addition to this feature, there are also small oscillations 
between these Andreev peaks. These oscillations are also results of the interference from all 
the interfaces. They have fixed period in terms of energy as demonstrated in Fig. 7(c). 

It is tempting to make an analogy between the NS superlattice and the diffraction grat- 
ings (each superconductor cell seems to be analogous to a slit in the grating). However, 
crucial differences do exist and make the interference in the NS superlattice system more 
complicated: (i) The Andreev reflection from each superconductor has to pass all other su- 
perconductors in front of it to get to the electrode where the interference happens, (ii) The 
incident wave that arrives at each superconductor has to pass all other superconductors in 
front of it. As a result, a precise dependence of the separation and the width on A and B 
are complicated and can only be computed numerically. Nevertheless, in a special limit in 
which A is large and B is very small, one may disregard B in computing the phase. The 
result is similar to that for the diffraction grating: Because the path between successive 
superconducting cells is 2AV , we have 



|C | 2 oc 



sin{N AVy' 2 



(33) 



sin(AV) 

Therefore, the positions of the main Andreev peaks are independent of the number of layers 
iVo and are located at ran /A with m being nonnegative integers. On the other hand, the 
minima are located at mir/NoA with m being nonnegative integers, and hence the width 
of the Andreev peak is tt/NqA. All of these results can be explicitly checked numerically. 
Fig. 7(d) is a simple demonstration in which the calculated dl/dV curve based on Eq . (|33"1) 
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is compared with that based on the exact calculation. One sees that Eq.(^) gives excellent 
results in the lower energy sector. Note that for large A, the Andreev peak gets sharper 
and becomes more like a zero-bias resonant conductance peak. This is very different from 
the other limit when A goes to zero, where the Andreev peak is a plateau, extending from 
zero-bias to V = A . 

We now analyze the IS superlattice. Essentially, the dl/dV curve for the IS superlattice 
has a similar behavior as shown in Fig. 7(e) in which an Andreev peak sits at the zero- 
bias, with reduced width. The difference is that the small oscillations now sit in a rising 
background. The reason for this rising behavior is due to that when V increases, the 
insulators becomes more transparent so that the conductance increases. The differential 
conductance shown Fig. 7(e) resembles a recent observed data in a underdoped BSCCO 
sample fllfl . This resemblance suggests that a IS superlattice is more appropriate to model 
the c-axis transport in the underdoped regime. 



IV. CONCLUSIONS AND ACKNOWLEDGMENTS 

In conclusion, we have investigated the effects of the superlattice structure on the single 
particle transport along the c-axis of the high Tc superconductors. Based on superlattice 
models that consist of metals/insulators and d-wave superconductors (NS/IS superlatives), 
we find that a crucial difference between the NS superlattice and the IS superlattice lies in 
the positions of the quasi-particle peaks in the measurement of density of state. In the large 
mass anisotropy limit of the metal for the NS superlattice, the density of state in the low 
energy section is still bulk d-wave like except that the position of the quasi-particle peak is 
reduced considerably, while for the IS superlattice, the quasi-particle peak remains at the gap 
value. The width of the Andreev peak at zero-bias in the planar junction measurement is also 
shown to be affected strongly by the superlattice structure. It is found that this width can 
be considerably reduced due to the destructive interference of the Andreev reflections from 
all the superconductors. In addition to this feature, there are also distinct oscillations with 
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smaller amplitudes extending out from the main Andreev peak. Most importantly, for the 
IS superlattice, these oscillations sit in a rising background, in consistent with experiments. 
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V. FIGURE CAPTIONS 

Fig.l The unit cell in a NS superlattice and the definition of angle <fi. We shall denote 
d = a + b and assume that the q x and q y directions are infinite. 

Fig. 2 (a) The one dimensional density of state for the case when 6 = 0, <fi = 0, m = m! 
and m^r = rris- Here A and B are chosen to be close the data for YBCO: A = 0.57, 
B = 0.21. It is seen that the density of state possess particle and hole symmetry, (b) 
The one dimensional density of state for the same parameters used in (a) except that now 
^s/ m N = 0.5. The asymmetry between particles and holes are evident. 

Fig. 3 (a) The three dimensional density of state for the homogeneous case when A = 
and B = 1. Here n(e) is measured by the unit k F /47r 2 d. (b) A similar plot for the case 
when A = 0.1,0.2 and B = 1. One sees that the quasi-particle peak moves in the subgap 
region, leaving a peak-like structure at e — 1. (c) A more realistic case when A = 0.57 and 
B = 0.21. 

Fig. 4 The three dimensional density of state in the infinite mass anisotropy limit 
[mjm! — > oo, A = 0.57, and B = 0.21) for the low energy section. The position of the 
quasi-particle is considerably smaller than Ao- 

Fig. 5 The three dimensional density of state for the SI superlattice in the simplest case 
when all the masses are the same and the other paramters are A = 0.57, and B = 0.21. 
Here V = 1.2ep and Bp = 10. One sees that the quasi-particle peak remains at e — 1. 

Fig. 6 Schematic plot of the setup for low contact resistance measurement. 

Fig. 7 (a) The calculated differential conductance for a NSN junction based on Eqs.(pO|) 
and (|3i~|) for different lengths of the superconductor. Here = 0, ef = 100, and dl/dV 
is measured in terms of 2N(EF)evpA. (b) The dl/dV curve (integrated over <f>) for a NS 
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superlattice with 10 layers. Here we assume m = m! = ms = m^, A = 0.57, B = 0.21, 
Z = 0, and ep = 20. (c) The energies of local maxima of the dl/dV versus their indices. 
Except a few steps in-between, the energy is linear in the index, indicating that the period 
of small oscillations in Fig.7(b) is a constant. The steps occurs when one jumps across the 
main Andreev peaks, (d) Solid line: The dl/dV curve for a NS superlattice with 15 layers 
for A = 1 and B = 0.01 . Here we assume that m = ml = ms = m^, Z = 0, = and 
Sf — 20. Dash line: A fitting based on Eq.([33|) gives excellent results in low energies, (e) 
The differential conductance for a SI superlattice with 10 layers. The parameters used here 
are A = 0.57, B = 0.21, Z = 0, 71 = 0.5, 72 = 2.87, e F = 20, V/A = 22. This curve is close 
to to the dl/dV curve in low temperature recently observed in Ref. JTTJ (see their Fig. 2(b)). 



21 




Fig. 1 



Fig. 2(a) 




Fig. 2(b) 

10.0 . r- — -n - ^ 



5.0 



0.0 



-7.0 



-5.0 



-3.0 



-1.0 



1.0 



3.0 



5.0 



7.0 



Fig.3(a) 




Fig. 3(b) 




Fig. 3(c) 




Fig. 4 




Fig. 5 




Fig. 6 



N 




V(x) 

A 



V 







N 



N 





N 



< d > 




V(l-((n-l)d+b)/ N d) 



V(l-n/ N n d) 




•N. 



-> x 



-a b 



nd (n+l)d 
n d+b 



N d 



Fig.7(a) 




Fig. 7(b) 




Fig.7(c) 



> 




120 



index of local maximum dl/dV 



Fig. 7(d) 



1.030 




0.990 

0.0 



2.0 

V 



4.0 



